#pragma once
#include <math.h>
#include "../params.h"

//Log of wave functions and derrivatives
inline double defaultPsi(double r)
{
	static const double curveFit = 0.63837943620287906921;

	double ri = 1.0 / (r - p4);
	double riPow2 = ri * ri;

	return -1.0 * ((p0 * riPow2 * ri + p1) * riPow2 + p2 * log(r)) + curveFit;
}

double defaultPsiPrime(double r)
{
  double rInv = 1.0 / (r - p4);
  double rInvPow3 = rInv * rInv * rInv;

  return (5.0 * p0 * rInvPow3 + 2.0 * p1) * rInvPow3 - p2 * rInv - p3;
}

double defaultPsiDoublePrime(double r)
{
  double rInv = 1.0 / (r - p4);
  double rInvPow2 = rInv * rInv;
  double rInvPow3 = rInvPow2 * rInv;

  return (-1.0 * (30.0 * p0 * rInvPow3 + 6.0 * p1) * rInvPow2 - p2) * rInvPow2;
}

//McMillan Wave functions and derrivatives
double mcMillanPsi(double r)
{
  return c00 * pow(r, c10) + c20 * r * r;
}

double mcMillanPsiPrime(double r)
{
  return c01 * pow(r, c11) + c21 * r;
}

double mcMillanPsiDoublePrime(double r)
{
  return c02 * pow(r, c12) + c22;
}

//Potential
double computePotential(double r)
{
	double x = r/rm;
	double xInverse = 1.0 / x;
	double xInversePow2 = xInverse * xInverse;
	double xInversePow6 = xInversePow2 * xInversePow2 * xInversePow2;
	double potential1 = A * exp(x * (beta * x - alphaPotential));
	double potential2 = (c6_pot + 
                       (c8_pot + c9_pot*xInverse + c10_pot * xInversePow2) * 
											xInversePow2) * xInversePow6;
	double potential3 = 1.0;

	if(x < D)
	{
		double tmp = (D / x) - 1.0;
		potential3 = exp(-1.0 * tmp * tmp);
	}

	return epsilon * (potential1 - (potential2 * potential3));
}
